Tema 7: Introducción a los métodos de análisis en R

Ir arriba

Métodos de análisis en R

El presente tema presenta algunos métodos de análisis en R.

Los temas cubiertos serán:

• Muestreo.
• Componentes principales.
• Análisis de Factores.
• Agrupamiento (cluster).
• Discriminante.
• Regresión.
• Regresión dicotómica.
• Series temporales.

Al finalizar la sesión usted tendrá una noción de una gama de análisis que se pueden realizar en R, esperando que incentive a la persona a seguir buscando otras formas, métodos y técnicas de análisis de la información.

1. Métodos de muestreo: selección de una muestra.

¿Qué es el muestreo?

El muestreo es el proceso de seleccionar un conjunto de individuos de una población con el fin de estudiarlos y poder caracterizar el total de la población. Dada la imposibilidad de estudiar o análisis a una población entera (censar), debemos pasar el una muestreo, o por el método de muestreo.

El muestreo se constituye de dos partes:

No sé cubrirá el tamaño de muestre (n), y sus distintos diseños (aleatorio, sistemático, estratificado, conglomerado o complejo, etc) pero veremos dos formas de seleccionar casos en R.

Muestra Aleatoria

Podemos tomar una muestra de tamaño n mediante la función “sample_n”. Tomemos una muestra de tamaño 50.

dim(iris)

## [1] 150   5

n <- iris %>% sample_n(50,replace=FALSE)
dim(n)

## [1] 50  5

Muestra según fracción de la población

También podemos seleccionar de forma aleatoria las observaciones pero de forma relativo o por una proporción.

Utilizamos la función “sample_frac”. Tomemos un 10% de la población.

dim(iris)

## [1] 150   5

n2 <- iris %>% sample_frac(0.1,replace=FALSE)
dim(n2)

## [1] 15  5

Otras funciones

R posee dos librerías las siguientes librerías para tratar los temas de muestreo:

• sampling
• survey

Mientras que sampling trabaja diseños aleatorios irrestrictos, sistemáticos y jerárquicos, survey se utiliza para diseños complejos o por conglomerados.

Tipo de análisis - NS.

Abordamos dos tipos de enfoque:

• Aprendizaje no supervisado y
• Aprendizaje supervisado.

Antes:

• ¿Qué es Machine Learning?
• ¿Cuál es la diferencia entre el Machine Learing y los métodos de Estadística?

Una definición de la web explica:

Aprendizaje no supervisado es un método donde un modelo se ajusta a las observaciones. Se distingue del Aprendizaje supervisado por el hecho de que no hay un conocimiento a priori. En el aprendizaje no supervisado, un conjunto de datos de objetos de entrada es tratado.

En mis palabras, diría:

• En el aprendizaje no supervisado, todas las variables tienen la misma relación entre ellas (relación simétrica). Ninguna busca explicar la relación o el comportamiento de las otras. Todas tienen el mismo nivel de explicación.
• En el aprendizaje supervisado una variable trata de ser explicada por las demás (relación asimétrica). Acá es donde se suele hablar de variable(s) dependiente(s) y variable(s) independiente(s).

Empecemos primero los análisis de aprendizaje no supervisado.

2. Análisis por Componentes Principales.

El análisis por componente principales hace parte de un grupo de análisis descriptivos multidimensionales llamados métodos factoriales.

Son métodos descriptivos, no se apoyan en modelos probabilísticos sino más bien de un modelo geométrico para mejorar la representación multidimencional.

A partir de una matriz rectangular de datos con p variables cuantitativas y n unidades, el análisis por componentes principales propone diversas representación geométricas para el entendimiento de los individuos y las variables.

Lo que se busca es ver si existe una estructura, no conocida a priori, para el conjunto de casos y variables, y así mejorar la interpretación.

Como todo método descriptivo, llevar a cabo un PCA no es un fin en sí. El PCA sirve para conocer mejor los datos, detectar valores sospechosos, y ayuda a formular hipótesis que se deben estudiar luego mediante modelos predictivos.

Para la aplicación del PCA en R, se pueden utilizar diversas librerías y funciones. Se exponen algunos ejemplos:

Veamos la estrucuta de las salidas del análisis por PCA.

Variancia explicada

Porcentaje de la variancia explicada

Veamos el porcentaje de variancia explicada por componente

Proporció de la varianza explicada

Proporció de la varianza explicada acumulada

Gráfico de sedimentación

Proyecciones

El PCA busca mediante proyección explicar a:

• Las observaciones o individuos.
• Las variables.
• La conjunción de observaciones-variables.

Individuos

Proyección de los individuos

Variables

Proyección de las variables

Individuos - Variables

Proyección de individuos y variables

3. Análisis de factores.

En el AF la atención se centra principalmente en las variables estadísticas y no en los individuos; es más un método de análisis multivariante que un método de análisis multidimensional.

El Análisis Factorial es un nombre genérico que se da a una clase de métodos estadísticos multivariantes cuyo propósito principal es definir la estructura subyacente en una matriz de datos. Generalmente hablando, aborda el problema de cómo analizar la estructura de las interrelaciones (correlaciones) entre un gran número de variables.

Al igual que el PCA, en el AF se busca reducir un número p a k variables. Sin embargo, el AF posee ciertas características que lo hacen el análisis ideal para conocer las estructura subyacente que se desea estudiar. Se dice que el AF es más elaborado, y se buscan conclusiones más contundentes, el PCA busca más aproximaciones a nivel gráfico del conjunto de datos. A esa reducción les llamamos Factores.

De forma matemática, el PCA es un AF pero más simple: sin extracción de vector propio, sin rotaciones, y sin transformaciones, además del gran supuesto practicamente irreal de la ortogonalidad de las variables o componentes principales en el PCA.

Desde un punto de vista práctico, el PCA es observacional, y el AF es un enfoque para un modelo explicito (estructura subyacente).

Mientras que el PCA busca el análisis de la variancia, el AF además busca el análisis de la covariancia y correlación entre las variables.

El objetivo típico en el análisis factorial es identificar las variables que están relacionadas entre sí, y separarlas de otras (una forma de agrupación de las variable). Los tipos de rotaciones ayudan a visualizar mejor este proceso.

Existen dos grandes tipos de modalidaes para el FA:

• Análisis Factorial Exploratorio.
• Análisis Factorial Confirmatorio.

Veamos unos usos del FA exploratorio.

Cantidad de Factores

Gráfico de sedimentación (Scree plot)

Opción 1

La forma más popular de conocer la cantidad óptima de factores es mediante el gráfico de sedimentación. Lo usual es ver dónde es que se lleva a cabo la caída más importante, o se dibuja un “codo”. A veces no es nada fácil el saber lo anterior. Personalmente considero que ante la poca claridad, hay que llevar a cabo un análisis interactivo del FA.

¿Cuántos factores?

Opción 2

R posee ciertos criterios que permiten “guiar” al analista en la elección del número de factor a elegir.

## Parallel analysis suggests that the number of factors =  3  and the number of components =  NA

##   noc naf nparallel nkaiser
## 1   5   1         5       5

Entonces, ¿cuántos componentes debemos elegir?

fa.1 <- fa(r=wine.1, nfactors = 3, rotate = "varimax")

res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)
fa.0 <- fa(r=personality, nfactors = 5, rotate = "varimax")
res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)
res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)
res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)

# Calculadon la  singularidad  y la comunalidad

loadings_distant = res1a$loadings[1,]
communality_distant = sum(loadings_distant^2); communality_distant

## [1] 0.4554525

uniqueness_distant = 1-communality_distant; uniqueness_distant

## [1] 0.5445475

Verificación de resultados

Contribuciones

Veamos el peso de cada variable a los Factores

Factor 1

Factor 2

Factor 3

Reporte del FA

Según la cantidad de factores, la función fa.diagram nos puede indicar cuál es la estructura factorial del análisis.

Factores en plano 2D

LLevado a cabo el proceso de elegir la rotación y luego el número de componentes, podemos verificar que desde un inicio hasta la solución del problema, la estructura factorial del FA posee una forma de explicar los facores en un plano 2D. Veamos esto

Inicio: todas las variables

Intermedio: aplicación del FA y el acomodo de las variables

Final: estructura final del FA

Creamos la estructura factorial mediante el análisis de las cargas factoriales

shy = rowMeans(cbind(personality$distant, personality$shy, personality$withdrw, personality$quiet))
outgoing = rowMeans(cbind(personality$talkatv, personality$outgoin, personality$sociabl))
hardworking = rowMeans(cbind(personality$hardwrk, personality$persevr, personality$discipl))
friendly = rowMeans(cbind(personality$friendl, personality$kind, personality$coopera, personality$agreebl, personality$approvn, personality$sociabl))
anxious = rowMeans(cbind(personality$tense, personality$anxious, personality$worryin))

# Combinando los factores y creando una nueva estructura de datos
combined_data = cbind(shy,outgoing,hardworking,friendly,anxious)
combined_data = as.data.frame(combined_data)

Puede que acá se estimen diversos modelos de FA, con diversas rotaciones, con tal de ver una mejorar separación, y hasta aglomeración de los factores.

4. Análisis de correspondencia.

Es una técnica descriptiva o exploratoria cuyo objetivo es resumir una gran cantidad de datos en un número reducido de dimensiones, con la menor pérdida de información posible. En esta línea, su objetivo es similar al de los métodos factoriales, salvo que en el caso del análisis de correspondencias el método se aplica sobre variables categóricas (nominales y ordinales).

El análisis de correspondencias simples se utiliza a menudo en la representación de datos que se pueden presentar en forma de tablas de contingencia de dos variables nominales u ordinales. Otras utilizaciones implican el tratamiento de tablas de proximidad o distancia entre elementos, y tablas de preferencias.

El análisis de correspondencias consiste en resumir la información presente en las filas y columnas de manera que pueda proyectarse sobre un subespacio reducido, y representarse simultáneamente los puntos fila y los puntos columna, pudiéndose obtener conclusiones sobre relaciones entre las dos variables nominales u ordinales de origen.

La extensión del análisis de correspondencias simples al caso de varias variables nominales (tablas de contingencia multidimensionales) se denomina Análisis de Correspondencias Múltiples, y utiliza los mismos principios generales que la técnica anterior. En general se orienta a casos en los cuales una variable representa ítems o individuos y el resto son variables cualitativas u ordinales que representan cualidades.

Librerías

Ejemplo de un CA simple

housetasks

##            Wife Alternating Husband Jointly
## Laundry     156          14       2       4
## Main_meal   124          20       5       4
## Dinner       77          11       7      13
## Breakfeast   82          36      15       7
## Tidying      53          11       1      57
## Dishes       32          24       4      53
## Shopping     33          23       9      55
## Official     12          46      23      15
## Driving      10          51      75       3
## Finances     13          13      21      66
## Insurance     8           1      53      77
## Repairs       0           3     160       2
## Holidays      0           1       6     153

res.ca <- CA(housetasks, graph = FALSE)
get_ca_row(res.ca)

## Correspondence Analysis - Results for rows
##  ===================================================
##   Name       Description                
## 1 "$coord"   "Coordinates for the rows" 
## 2 "$cos2"    "Cos2 for the rows"        
## 3 "$contrib" "contributions of the rows"
## 4 "$inertia" "Inertia of the rows"

get_ca_col(res.ca)

## Correspondence Analysis - Results for columns
##  ===================================================
##   Name       Description                   
## 1 "$coord"   "Coordinates for the columns" 
## 2 "$cos2"    "Cos2 for the columns"        
## 3 "$contrib" "contributions of the columns"
## 4 "$inertia" "Inertia of the columns"

#--- para las filas:
fviz_contrib(res.ca, choice = "row", axes = 1)

#--- para las columnas:
fviz_contrib(res.ca, choice = "col", axes = 1)

#--- para las variables en las filas:
fviz_ca_row(res.ca, repel = TRUE)

#--- para las variables en las columnas:
fviz_ca_col(res.ca)

#--- y la visualización conjunta:
fviz_ca_biplot(res.ca, repel = TRUE)

5. Análisis de agrupación o Cluster Analysis.

El Análisis por Conglomerados (Cluster Analysis), es una técnica estadística multivariante que agrupa elementos (sujetos) tratando de lograr la máxima homogeneidad en cada grupo y la mayor diferencia entre los grupos.

Es un método basado en criterios geométricos y se utiliza fundamentalmente como una técnica exploratoria, descriptiva pero no explicativa. Se fundamenta en teorías de optimizaciones numéricas, y no en optimizaciones algebraicas.

Las soluciones no son únicas, en la medida en que la pertenencia al conglomerado para cualquier número de soluciones depende de muchos elementos del procedimiento elegido.

La solución del Cluster depende totalmente de las variables utilizadas, la adición o destrucción de variables relevantes puede tener un impacto substancial sobre la solución resultante.

Los algoritmos de formación de conglomerados se agrupan en diversas categorías, y existen muchas modalidades. El presente capítulo estudia los modalidades de clusters jerárquicos y los de k-medias.

El objetivo del análisis de cluster es, a partir de “𝑛” individuos (𝑛 grupos iniciales) buscar agrupar el conjunto de individios en dos o más grupos basados en la similitud de esos objetos según una serie especificada de características (variables).

Finalmente, siempre hay que tomar en cuenta los mismos tres aspectos que podrían influir fuertemente en el análisis de los datos:

• Valores extremos
• Medición de la similitud
• Estandarización de los datos

Veremos algunos ejemplos de los algoritmos presentados.

Funciones y librerías del cluster analysis

En el entorno de programación R existen múltiples paquetes que implementan algoritmos de clustering y funciones para visualizar sus resultados. En este documento se emplean los siguientes:

stats: contiene las funciones:

• dist() para calcular matrices de distancias
• kmeans()
• hclust()
  
• cuttree() para crear los clusters y
• plot.hclust() para visualizar los resultados.
• luster
• mclust: contienen múltiples algoritmos de clustering y métricas para evaluarlos.
• factoextra: extensión basada en ggplot2 para crear visualizaciones de los resultados de clustering y su evaluación.
• dendextend: extensión para la customización de dendrogramas.

Para evitar problemas de la magnitud de las variables, recordar SIEMPRE escalar las variables…

Veamos la estructura de algunas de funciones anteriores:

scale

La estructura es :

#scale(x, center = TRUE, scale = TRUE)

kmeans

La estructura es :

#kMeans(x, centers, iter.max=10, num.seeds=10)

fviz_cluster

La estructura es :

# fviz_cluster(object, data = NULL, choose.vars = NULL, stand = TRUE,
#  axes = c(1, 2), geom = c("point", "text"), repel = FALSE,
#  show.clust.cent = TRUE, ellipse = TRUE, ellipse.type = "convex",
#  ellipse.level = 0.95, ellipse.alpha = 0.2, shape = NULL,
#  pointsize = 1.5, labelsize = 12, main = "Cluster plot", xlab = NULL,
#  ylab = NULL, outlier.color = "black", outlier.shape = 19,
#  ggtheme = theme_grey(), ...)

hclust

La estructura es :

# hclust(d, method = "complete", members = NULL)
# S3 method for hclust
# plot(x, labels = NULL, hang = 0.1, check = TRUE,
#      axes = TRUE, frame.plot = FALSE, ann = TRUE,
#      main = "Cluster Dendrogram",
#      sub = NULL, xlab = NULL, ylab = "Height", …)

Cluster por k medias

El set de datos USArrests contiene información sobre el número de delitos (asaltos, asesinatos y secuestros) junto con el porcentaje de población urbana para cada uno de los 50 estados de USA. Se pretende estudiar si existe una agrupación subyacente de los estados empleando K-means-clustering.

El paquete factoextra creado contiene funciones que facilitan en gran medida la visualización y evaluación de los resultados de clustering. Si se emplea K-means-clustering con distancia euclídea hay que asegurarse de que las variables empleadas son de tipo continuo, ya que trabaja con la media de cada una de ellas.

Veamos una representación de los datos, para 4 clusters (K=4)

datos <- scale(USArrests)

set.seed(123)
km_clusters <- kmeans(x = datos, centers = 4, nstart = 50)

# Las funciones del paquete factoextra emplean el nombre de las filas del
# dataframe que contiene los datos como identificador de las observaciones.
# Esto permite añadir labels a los gráficos.
fviz_cluster(object = km_clusters, data = datos, show.clust.cent = TRUE,
             ellipse.type = "euclid", star.plot = TRUE, repel = TRUE) +
  labs(title = "Resultados clustering K-means") +
  theme_bw() +
  theme(legend.position = "none")

Clusters jerárquicos

Hierarchical clustering es una alternativa a los métodos de partitioning clustering que no requiere que se pre-especifique el número de clusters. Los métodos que engloba el hierarchical clustering se subdividen en dos tipos dependiendo de la estrategia seguida para crear los grupos:

Agglomerative clustering (bottom-up): el agrupamiento se inicia en la base del árbol, donde cada observación forma un cluster individual. Los clusters se van combinado a medida que la estructura crece hasta converger en una única “rama” central.

Divisive clustering (top-down): es la estrategia opuesta al agglomerative clustering, se inicia con todas las observaciones contenidas en un mismo cluster y se suceden divisiones hasta que cada observación forma un cluster individual.

En ambos casos, los resultados pueden representarse de forma muy intuitiva en una estructura de árbol llamada dendrograma.

En este se debe antes transformar los datos a distancias

# Crear las distancia y definir el método 
dd <- dist(scale(USArrests), method = "euclidean")
hc <- hclust(dd, method = "ward.D2")

Veamos diversas modalidaes

Tipos de clusters jerárquicos

Referencia: http://www.sthda.com/english/wiki/beautiful-dendrogram-visualizations-in-r-5-must-known-methods-unsupervised-machine-learning

Mediante el plot.hclust(), el más sencillo

plot(hc)

# Put the labels at the same height: hang = -1
plot(hc, hang = -1, cex = 0.6)

Mediante el plot.dendogram(),

# Convert hclust into a dendrogram and plot
hcd <- as.dendrogram(hc)
# Default plot
plot(hcd, type = "rectangle", ylab = "Height")

# Zoom in to the first dendrogram
plot(hcd, xlim = c(1, 20), ylim = c(1,8))

# Triangle plot
plot(hcd, type = "triangle", ylab = "Height")

Mediante el Phylogenetic trees

Este presenta diversas formas de representar los clusters

# Default plot
plot(as.phylo(hc), cex = 0.6, label.offset = 0.5)

# Cladogram
plot(as.phylo(hc), type = "cladogram", cex = 0.6, 
     label.offset = 0.5)

# Unrooted
plot(as.phylo(hc), type = "unrooted", cex = 0.6,
     no.margin = TRUE)

# Fan
plot(as.phylo(hc), type = "fan")

# Radial
plot(as.phylo(hc), type = "radial")

# Change the appearance
# change edge and label (tip)
plot(as.phylo(hc), type = "cladogram", cex = 0.6,
     edge.color = "steelblue", edge.width = 2, edge.lty = 2,
     tip.color = "steelblue")

Tipo de análisis - S.

Ahora pasemos a los análisis supervisados.

En el aprendizaje supervisado una variable trata de ser explicada por las demás (relación asimétrica, variables dependientes e independientes).

Clasificación / análisis discriminante.

La clasificación se centra en identificar a cuál de un conjunto de categorías (subpoblaciones) pertenece una nueva observación, sobre la base de un conjunto de datos de formación que contiene observaciones (o instancias) cuya categoría de miembros es conocida.

En el análisis de aprendizaje automático, para datos transversales, existen dos grandes categorias: la clasificación y la predicción.

Por ahora veremos la clasificación, luego en la regresión veremos la predicción, y volveremos a ver la clasifiación cuando vemos la regresión dicotómica.

6. Análisis discriminante

El análisis discriminante es una técnica multivariante cuya finalidad es describir las diferencias significativas entre 𝑘 grupos de objetos (𝑘 > 1) sobre los que se observan (variable de clasificación).

En caso de que estas diferencias existan, intentará explicar en qué sentido se dan y proporcionar procedimientos de asignación sistemática de nuevas observaciones con grupo desconocido a uno de los grupos analizados, utilizando para ello sus valores en las𝑘variables clasificadoras.

Podemos ver este procedimiento como un modelo de predicción de una variable respuesta categórica (variable grupo) a partir de 𝑘 variables explicativas generalmente continuas (variables clasificatorias).

El Análisis Discriminante consiste en buscar nuevas variables (variables discriminantes), que separan lo mejor posible los grupos de proyección en las 𝑘 variables, de las 𝑛 observaciones, para poder predecir nuevas observaciones i.

En la clasificación se llevan a cado dos procesos fundamentales: construir la regla o ecuación de clasificación, y verificar la calidad del proceso, para luego introducir nuevos casos..

En la construcción de la regla de decisión, veremos un forma lineal y otro cuadrática.

Librerías y funciones

Utilizaremos dos funciones lda() y qda() de la librería MASS.

lda() (MASS)

La estructura es :

 lda(x, …)
  S3 method for formula
 lda(formula, data, …, subset, na.action)

 S3 method for default
 lda(x, grouping, prior = proportions, tol = 1.0e-4,
    method, CV = FALSE, nu, …)

 S3 method for data.frame
 lda(x, …)

  S3 method for matrix
  lda(x, grouping, …, subset, na.action)

Ver el enlace: https://www.rdocumentation.org/packages/MASS/versions/7.3-51.1/topics/lda

qda()

La estructura es :

qda(x, …)
S3 method for formula
qda(formula, data, …, subset, na.action)

S3 method for default
qda(x, grouping, prior = proportions,
   method, CV = FALSE, nu, …)

S3 method for data.frame
qda(x, …)

S3 method for matrix
qda(x, grouping, …, subset, na.action)

Ver el enlace: https://www.rdocumentation.org/packages/MASS/versions/7.3-51.1/topics/qda

Partición del archivo de datos

Pera evaluar la regla de decisión, se suele particionar el archivo de datos en entramiento y validación. Esto es un procedimiento típico de los análisis supervisados, y también en la minería de datos.

training_sample <- sample(c(TRUE, FALSE), nrow(iris), replace = T, prob = c(0.6,0.4))

train <- iris[training_sample, ]
test <- iris[!training_sample, ]

Dimensiones del set de entrenamiento

## [1] 89  5

Dimensiones del set de validación

## [1] 61  5

Estimación del lda y qda

Estimación del lda

lda.iris <- lda(Species ~ ., train)
lda.iris

## Call:
## lda(Species ~ ., data = train)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3258427  0.3483146  0.3258427 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa         5.051724    3.403448     1.462069   0.2344828
## versicolor     5.974194    2.783871     4.290323   1.3709677
## virginica      6.541379    2.979310     5.517241   1.9965517
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8752702 -0.72840217
## Sepal.Width   1.4177443  3.31503628
## Petal.Length -2.0953556  0.07481247
## Petal.Width  -3.0258487  1.42440523
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9911 0.0089

names(lda.iris)

##  [1] "prior"   "counts"  "means"   "scaling" "lev"     "svd"     "N"      
##  [8] "call"    "terms"   "xlevels"

Enfoquemos en la salida, la parte de los coeficientes (podemos verlo vediante el “scaling”)

¿Cuál sería entonces las ecuaciones de ambas reglas?

LD1 ? LD2 ?

Estimación del qda

qda.iris <- qda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, train)
qda.iris

## Call:
## qda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
##     data = train)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3258427  0.3483146  0.3258427 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa         5.051724    3.403448     1.462069   0.2344828
## versicolor     5.974194    2.783871     4.290323   1.3709677
## virginica      6.541379    2.979310     5.517241   1.9965517

Analisando la función qda, no se programó la ecuación del qda…

Análisis de Regla de discriminación

Una forma de ver la efectividad de las reglas de separarción en poner en un plano 2D los casos con sus respectivas reglas si hubiese 2.

plot(lda.iris, col = as.integer(train$Species))

También se puede corroborar mediante el análisis de la posición y variabilidad por categoría de clasificación, en cada regla. El caso de la RD1:

lda.iris <- lda(Species ~ ., train)
lda.iris

## Call:
## lda(Species ~ ., data = train)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3258427  0.3483146  0.3258427 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa         5.051724    3.403448     1.462069   0.2344828
## versicolor     5.974194    2.783871     4.290323   1.3709677
## virginica      6.541379    2.979310     5.517241   1.9965517
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8752702 -0.72840217
## Sepal.Width   1.4177443  3.31503628
## Petal.Length -2.0953556  0.07481247
## Petal.Width  -3.0258487  1.42440523
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9911 0.0089

plot(lda.iris, dimen = 1, type = "b")

Tengo un error para la 2nd dimensión… creo que debe ser porque explica muy poco.

#plot(lda.iris, dimen = 2, type = "both")

Otra forma de ver la efectividad de las reglas es mediante la proyección de la regla y sus inviduos para ciertas variables.

Veamos que tanto separa o discrimina la primera regla

Y para la 2nda regla:

¿Qué podemos concluir?

Los gráficos de partición y las proyecciones

Particiones clásiscas con partimat

El uso de la función partimat del paquete klaR proporciona una forma alternativa de trazar las funciones discriminantes lineales. partimat emite una serie de gráficos para cada combinación de dos variables. Piense en cada gráfico como una vista diferente de los mismos datos. Las regiones coloreadas delinean cada área de clasificación. Se predice que cualquier observación que caiga dentro de una región será de una clase específica. Cada gráfico también incluye la tasa de error aparente para esa vista de los datos.

Partición lineal

Partición cuadrática

Particiones partimat cambiando el color

La función partimat() del paquete klar permite representar los límites de clasificación de un modelo discriminante lineal o cuadrático para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.

Partición lineal

Partición cuadrática

Proyección del modelo: variables y predicciones

Podemos verificar la calidad del análisis mediante las proyeccioes de las variables y las predicciones.

Variables: podemos ver los resultados del modelo de la separación de los casos en los variables.

¿Qué podemos decir?

Predicciones.

Podemos para los valores de las reglas verificar si los individuos se están o no separando o aglomerando en cierto parte.

LD1: ¿qué podemos decir?

LD2: ¿qué podemos decir?

Las tablas de confusión

Set de entrenamiento

LDA

##             
##              setosa versicolor virginica
##   setosa         29          0         0
##   versicolor      0         30         1
##   virginica       0          1        28

–> Recomendable tenerlo mejor en %

¿Qué podemos concluir del modelo para el entrenamiento?

QDA

##             
##              setosa versicolor virginica
##   setosa         29          0         0
##   versicolor      0         30         1
##   virginica       0          1        28

–> Recomendable tenerlo mejor en %

¿Qué podemos concluir del modelo para el entrenamiento? Y comparando LDA y QDA

Set de validación

LDA

##             
##              setosa versicolor virginica
##   setosa         21          0         0
##   versicolor      0         18         0
##   virginica       0          1        21

–> Recomendable tenerlo mejor en %

¿Qué podemos concluir del modelo para el entrenamiento?

QDA

##             
##              setosa versicolor virginica
##   setosa         21          0         0
##   versicolor      0         18         0
##   virginica       0          1        21

–> Recomendable tenerlo mejor en %

¿Qué podemos concluir del modelo para el entrenamiento? Y comparando LDA y QDA

7. Regresión.

El objetivo es mostrar los principales comandos en R para generar:

• Regresió lineal simple
• Regresión lineal múltiple.

De igual forma se verán ciertas etapas de análisis que componen a estas variantes de la regresión.

Datos

El dataset Boston del paquete MASS recoge la mediana del valor de la vivienda en 506 áreas residenciales de Boston. Junto con el precio, se han registrado 13 variables adicionales.

tail(Boston)

##        crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
## 501 0.22438  0  9.69    0 0.585 6.027 79.7 2.4982   6 391    19.2 396.90 14.33
## 502 0.06263  0 11.93    0 0.573 6.593 69.1 2.4786   1 273    21.0 391.99  9.67
## 503 0.04527  0 11.93    0 0.573 6.120 76.7 2.2875   1 273    21.0 396.90  9.08
## 504 0.06076  0 11.93    0 0.573 6.976 91.0 2.1675   1 273    21.0 396.90  5.64
## 505 0.10959  0 11.93    0 0.573 6.794 89.3 2.3889   1 273    21.0 393.45  6.48
## 506 0.04741  0 11.93    0 0.573 6.030 80.8 2.5050   1 273    21.0 396.90  7.88
##     medv
## 501 16.8
## 502 22.4
## 503 20.6
## 504 23.9
## 505 22.0
## 506 11.9

• crim: ratio de criminalidad per cápita de cada ciudad.
• zn: Proporción de zonas residenciales con edificaciones de más de 25.000 pies cuadrados.
• indus: proporción de zona industrializada.
• chas: Si hay río en la ciudad (= 1 si hay río; 0 no hay).
• nox: Concentración de óxidos de nitrógeno (partes per 10 millón).
• rm: promedio de habitaciones por vivienda.
• age: Proporción de viviendas ocupadas por el propietario construidas antes de 1940.
• dis: Media ponderada de la distancias a cinco centros de empleo de Boston.
• rad: Índice de accesibilidad a las autopistas radiales.
• tax: Tasa de impuesto a la propiedad en unidades de $10,000.
• ptratio: ratio de alumnos/profesor por ciudad.
• black: 1000(Bk - 0.63)^2 donde Bk es la proporción de gente de color por ciudad.
• lstat: porcentaje de población en condición de pobreza.
• medv: Valor mediano de las casas ocupadas por el dueño en unidades de $1000s.

Regresión lineal simple

Análisis descriptivos de los datos

summary(Boston)

##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

Veamos las distribuciones por variable:

multi.hist(x = Boston[,1:3], dcol = c("blue","red"), dlty = c("dotted", "solid"),
           main = "")

multi.hist(x = Boston[,5:9], dcol = c("blue","red"), dlty = c("dotted", "solid"),
           main = "")

multi.hist(x = Boston[,10:14], dcol = c("blue","red"),
           dlty = c("dotted", "solid"), main = "")

vemos que los datos no poseen total normalidad en las variables.

Regresión lineal simple

Se pretende predecir el valor de la vivienda en función del porcentaje de pobreza de la población. Empleando la función lm() se genera un modelo de regresión lineal por mínimos cuadrados en el que la variable respuesta es medv y el predictor lstat.

modelo_simple <- lm(data = Boston,formula = medv ~ lstat)
modelo_simple

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Coefficients:
## (Intercept)        lstat  
##       34.55        -0.95

names(modelo_simple)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

summary(modelo_simple)

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.168  -3.990  -1.318   2.034  24.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.55384    0.56263   61.41   <2e-16 ***
## lstat       -0.95005    0.03873  -24.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared:  0.5441, Adjusted R-squared:  0.5432 
## F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

En la información devuelta por el summary se observa que el p-value del estadístico F es muy pequeño, indicando que al menos uno de los predictores del modelo está significativamente relacionado con la variable respuesta. ¿Cómo sería en palabras del presente ejemplo?

La creación de un modelo de regresión lineal simple suele acompañarse de una representación gráfica superponiendo las observaciones con el modelo. Además de ayudar a la interpretación, es el primer paso para identificar posibles violaciones de las condiciones de la regresión lineal.

attach(Boston)
plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
abline(modelo_simple, lwd = 3, col = "red")

La representación gráfica de las observaciones muestra que la relación entre ambas variables estudiadas no es del todo lineal, lo que apunta a que otro tipo de modelo podría explicar mejor la relación. Aun así la aproximación no es mala.

Una vez generado el modelo, es posible predecir el valor de la vivienda sabiendo el estatus de la población en la que se encuentra. Toda predicción tiene asociado un error y por lo tanto un intervalo. Es importante diferenciar entre dos tipos de intervalo:

Intervalo de confianza: Devuelve un intervalo para el valor promedio de todas las viviendas que se encuentren en una población con un determinado porcentaje de pobreza, supóngase lstat=10.

predict(object = modelo_simple, newdata = data.frame(lstat = c(10)),
        interval = "confidence", level = 0.95)

##        fit      lwr      upr
## 1 25.05335 24.47413 25.63256

Intervalo de predicción: Devuelve un intervalo para el valor esperado de una vivienda en particular que se encuentre en una población con un determinado porcentaje de pobreza.

predict(object = modelo_simple, newdata = data.frame(lstat = c(10)),
        interval = "prediction", level = 0.95)

##        fit      lwr      upr
## 1 25.05335 12.82763 37.27907

Una de las mejores formas de confirmar que las condiciones necesarias para un modelo de regresión lineal simple por mínimos cuadrados se cumplen es mediante el estudio de los residuos del modelo.

En R, los residuos se almacenan dentro del modelo bajo el nombre de residuals. R genera automáticamente los gráficos más típicos para la evaluación de los residuos de un modelo.

par(mfrow = c(1,2))
plot(modelo_simple)

par(mfrow = c(1,1))

Otra forma de identificar las observaciones que puedan ser outliers o puntos con alta influencia (leverage) es emplear las funciones rstudent() y hatvalues().

plot(x = modelo_simple$fitted.values, y = abs(rstudent(modelo_simple)),
     main = "Absolute studentized residuals vs predicted values", pch = 20,
     col = "grey30")
abline(h = 3, col = "red")

plot(hatvalues(modelo_simple), main = "Medición de leverage", pch = 20)
# Se añade una línea en el threshold de influencia acorde a la regla
# 2.5x((p+1)/n)
abline(h = 2.5*((dim(modelo_simple$model)[2]-1 + 1)/dim(modelo_simple$model)[1]),
       col = "red")

En este caso muchos de los valores parecen posibles outliers o puntos con alta influencia porque los datos realmente no se distribuyen de forma lineal en los extremos. Deberíamos evaluar otros modelos y otras variables.

Regresión múltiple

Se desea generar un modelo que permita explicar el precio de la vivienda de una población empleando para ello cualquiera de las variables disponibles en el dataset Boston y que resulten útiles en el modelo.

R permite crear un modelo con todas las variables incluidas en un data.frame de la siguiente forma:

modelo_multiple <- lm(formula = medv ~ ., data = Boston)
# También se pueden especificar una a una 
summary(modelo_multiple)

## 
## Call:
## lm(formula = medv ~ ., data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.595  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
## crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
## zn           4.642e-02  1.373e-02   3.382 0.000778 ***
## indus        2.056e-02  6.150e-02   0.334 0.738288    
## chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
## nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
## rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
## age          6.922e-04  1.321e-02   0.052 0.958229    
## dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
## rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
## tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
## ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
## black        9.312e-03  2.686e-03   3.467 0.000573 ***
## lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
## F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

En el summary se puede observar que algunos predictores tienen p-values muy altos, sugiriendo que no contribuyen al modelo por lo que deben ser excluidos, por ejemplo age e indus. La exclusión de predictores basándose en p-values no es aconsejable, en su lugar se recomienda emplear métodos de best subset selection, stepwise selection (forward, backward e hybrid) o Shrinkage/regularization.

step(modelo_multiple, direction = "both", trace = 0)

## 
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat, data = Boston)
## 
## Coefficients:
## (Intercept)         crim           zn         chas          nox           rm  
##   36.341145    -0.108413     0.045845     2.718716   -17.376023     3.801579  
##         dis          rad          tax      ptratio        black        lstat  
##   -1.492711     0.299608    -0.011778    -0.946525     0.009291    -0.522553

La selección de predictores empleando stepwise selection (hybrid/doble) ha identificado como mejor modelo el formado por los predictores crim, zn, chas, nox, rm, dis, rad, tax, ptratio, black, lstat.

modelo_multiple <- lm(formula = medv ~ crim + zn + chas +  nox + rm +  dis +
                      rad + tax + ptratio + black + lstat, data = Boston)
# También se pueden indicar todas las variables de un data.frame y exluir algunas
# modelo_multiple <- lm(formula = medv~. -age -indus, data = Boston)
summary(modelo_multiple)

## 
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5984  -2.7386  -0.5046   1.7273  26.2373 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
## crim         -0.108413   0.032779  -3.307 0.001010 ** 
## zn            0.045845   0.013523   3.390 0.000754 ***
## chas          2.718716   0.854240   3.183 0.001551 ** 
## nox         -17.376023   3.535243  -4.915 1.21e-06 ***
## rm            3.801579   0.406316   9.356  < 2e-16 ***
## dis          -1.492711   0.185731  -8.037 6.84e-15 ***
## rad           0.299608   0.063402   4.726 3.00e-06 ***
## tax          -0.011778   0.003372  -3.493 0.000521 ***
## ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
## black         0.009291   0.002674   3.475 0.000557 ***
## lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
## F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

Supuestos

En los modelos de regresión lineal con múltiples predictores, además del estudio de los residuos vistos en el modelo simple, es necesario descartar colinealidad o multicolinealidad entre variables.

par(mfrow = c(1,2))
plot(modelo_multiple)

par(mfrow = c(1,1))

Para la colinealidad se recomienda calcular el coeficiente de correlación entre cada par de predictores incluidos en el modelo:

require(corrplot)
corrplot.mixed(corr = cor(Boston[,c("crim", "zn", "nox", "rm", "dis", "rad", 
                                    "tax", "ptratio", "black", "lstat", "medv")],
                          method = "pearson"))

El análisis muestra correlaciones muy altas entre los predictores rad y tax (positiva) y entre dis y nox (negativa).

attach(Boston)

## The following objects are masked from Boston (pos = 3):
## 
##     age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
##     rm, tax, zn

par(mfrow = c(2,2))
plot(x = tax, y = rad, pch = 20)
plot(x = tax, y = nox, pch = 20)
plot(x = dis, y = nox, pch = 20)
plot(x = medv, y = rm, pch = 20)

par(mfrow = c(1,1))

Si la correlación es alta y por lo tanto las variables aportan información redundante, es recomendable analizar si el modelo mejora o no empeora excluyendo alguno de estos predictores.

Para el estudio de la multicolinealidad una de las medidas más utilizadas es el factor de inflación de varianza VIF. Puede calcularse mediante la función vif() del paquete car.

require(car)

## Loading required package: car

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:boot':
## 
##     logit

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

vif(modelo_multiple)

##     crim       zn     chas      nox       rm      dis      rad      tax 
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386 
##  ptratio    black    lstat 
## 1.757681 1.341559 2.581984

Los indices VIF son bajos o moderados, valores entre 5 y 10 indican posibles problemas y valores mayores o iguales a 10 se consideran muy problemáticos.

De igual forma, podemos realizar la predicción para el modelo múltiple pero indicando el valor a predir para cada variable.

Regresión Polinomial: incorporar no-linealidad a los modelos lineales

La Regresión Polinomial, aunque permite describir relaciones no lineales, se trata de un modelo lineal en el que se incorporan nuevos predictores elevando el valor de los ya existentes a diferentes potencias.

Cuando se intenta predecir el valor de la vivienda en función del estatus de la población, el modelo lineal generado no se ajusta del todo bien debido a que las observaciones muestran una relación entre ambas variables con cierta curvatura.

attach(Boston)

## The following objects are masked from Boston (pos = 5):
## 
##     age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
##     rm, tax, zn

## The following objects are masked from Boston (pos = 6):
## 
##     age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
##     rm, tax, zn

plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
abline(modelo_simple, lwd = 3, col = "red")

La curvatura descrita apunta a una posible relación cuadrática, por lo que un polinomio de segundo grado podría capturar mejor la relación entre las variables. En R se pueden generar modelos de regresión polinómica de diferentes formas:

Identificando cada elemento del polinomio: modelo_pol2 <- lm(formula = medv ~ lstat + I(lstat^2), data = Boston) El uso de I() es necesario ya que el símbolo ^ tiene otra función dentro de las formula de R.

Con la función poly(): lm(formula = medv ~ poly(lstat, 2), data = Boston)

modelo_pol2 <- lm(formula = medv ~ poly(lstat, 2), data = Boston)
summary(modelo_pol2)

## 
## Call:
## lm(formula = medv ~ poly(lstat, 2), data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2834  -3.8313  -0.5295   2.3095  25.4148 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.5328     0.2456   91.76   <2e-16 ***
## poly(lstat, 2)1 -152.4595     5.5237  -27.60   <2e-16 ***
## poly(lstat, 2)2   64.2272     5.5237   11.63   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.524 on 503 degrees of freedom
## Multiple R-squared:  0.6407, Adjusted R-squared:  0.6393 
## F-statistic: 448.5 on 2 and 503 DF,  p-value: < 2.2e-16

El p-value próximo a 0 del predictor cuadrático de lstat indica que contribuye a mejorar el modelo.

attach(Boston)

## The following objects are masked from Boston (pos = 3):
## 
##     age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
##     rm, tax, zn

## The following objects are masked from Boston (pos = 6):
## 
##     age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
##     rm, tax, zn

## The following objects are masked from Boston (pos = 7):
## 
##     age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
##     rm, tax, zn

plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
points(lstat, fitted(modelo_pol2), col = 'red', pch = 20)

A la hora de comparar dos modelos se pueden evaluar sus R2. En este caso el modelo cuadrático es capaz de explicar un 64% de variabilidad frente al 54% del modelo lineal.

Dado que un polinomio de orden n siempre va a estar anidado a uno de orden n+1, se pueden comparar modelos polinómicos dentro un rango de grados haciendo comparaciones secuenciales.

anova(modelo_simple, modelo_pol2)

## Analysis of Variance Table
## 
## Model 1: medv ~ lstat
## Model 2: medv ~ poly(lstat, 2)
##   Res.Df   RSS Df Sum of Sq     F    Pr(>F)    
## 1    504 19472                                 
## 2    503 15347  1    4125.1 135.2 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p-value obtenido para el estadístico F confirma que el modelo cuadrático es superior.

8. Regresión dicotómica.

Si en vez de una variable continua tenemos una variable categórica, y esa desea ser explicada por un conjunto de predictores o variables independienntes, podemos aplicar una regresión dicotómica.

Existen dos grandes tipos:

• Logística
• Probit

Datos

• Descripción: registros con información de 79 muestras de orina con características físicas.
• Variables: están registradas las siguientes variables:
  • r: indicador de presencia de oxalato de calcio. (Variable Objetivo)
    • 0: ausencia de oxalato en orina.
    • 1: presencia de oxalato en orina.
  • gravity: gravedad específica de la orina.
  • ph: pH de la orina.
  • osmo: osmolaridad de la orina.
  • cond: conductivicad de la orina.
  • urea: concentración de urea en la orina.
  • calc: concentración de calcio (milimoles por litro).
• Problema: determinar si una o más características físicas de la orina están relacionadas con la presencia de oxalatos. Además, generar un modelo capaz de clasificar pacientes con presencia de cristales que podrían ser causantes de patologías (cálculos renales).

datos <- read.csv("C:/Users/oscar/Desktop/R --- SAF/Tema 7/Orina.csv")

# Conversión de variable objetivo a factor
datos$r <- as.factor(datos$r)

# Imprimiendo datos
tail(datos)

##    r gravity   ph osmo cond urea  calc
## 74 1   1.022 5.09  736 19.8  418  8.53
## 75 1   1.025 7.90  721 23.6  301  9.04
## 76 1   1.017 4.81  410 13.3  195  0.58
## 77 1   1.024 5.40  803 21.8  394  7.82
## 78 1   1.016 6.81  594 21.4  255 12.20
## 79 1   1.015 6.03  416 12.8  178  9.39

La regresión logística

• La regresión logística analiza datos con distribución binomial de la forma:

\[Y_i \sim\ B(p_i,\ n_i),\ para\ i=1,...,m\]

• En la expresión anterior \(p_i\) hace referencia a la probabilidad de éxito (probabilidad de que ocurra el evento bajo estudio) y \(n_i\) determina el número de ensayos tipo Bernoulli. El número de ensayos es conocido, sin embargo, la probabilidad del éxito se desconoce.
• Se debe cumplir que la respuesta esté acotada entre 0 y 1, es decir, que el resultado siempre será positivo, además de ser inferior a 1.
• El exponencial (\(e\)) de cualquier valor (\(x\)) es siempre positivo y, cualquier número divivido entre la cantidad más uno (\(x+1\)) siempre será menor que 1. Bajo estas dos premisas se puede expresar la siguiente probabilidad condicional (función logística):

\[p(Y =1\ |\ X)=\frac{e^{(\beta_0+\beta_1x)}}{e^{(\beta_0+\beta_1x)}+1}\]

• Para facilitar el cálculo escribimos \(p(Y =1\ |\ X)\) como \(p(X)\):

\[p(X) = \frac{e^{(\beta_0+\beta_1x)}}{e^{(\beta_0+\beta_1x)}+1}\\ p(e^{(\beta_0+\beta_1x)}+1) = e^{(\beta_0+\beta_1x)}\\ p \times e^{(\beta_0+\beta_1x)}\ +\ p = e^{(\beta_0+\beta_1x)}\\ p = e^{(\beta_0+\beta_1x)}\ -\ p \times e^{(\beta_0+\beta_1x)}\\ p = e^{(\beta_0+\beta_1x)}(1-p)\\ \frac{p}{1-p} = e^{(\beta_0+\beta_1x)}\]

• Los logits (función de enlace) de las probabilidades binomiales desconocidas, es decir, los logaritmos de la razón de momios (odds ratio) son modelados como una función lineal de los \(X_i\):

\[ln(\frac{p}{1-p}) = \beta_0+\beta_1x\]

• Esta función de enlace es conocida como sigmoide y limita su rango de probabilidades entre 0 y 1.

La regresión probit

• La regresión probit permite analizar datos con respuesta ordinal o con distribución binomial (respuestas dicotómicas) de la forma:

\[Y_i \sim\ B(p_i,\ n_i),\ para\ i=1,...,m\]

• El marco conceptual del modelo probit puede ser expresado de la siguiente manera:

\[p(Y = 1|X)=\ \Phi(X^{T}\beta)\]

• Donde \(p(Y = 1|X)\) denota la probabilidad, \(\Phi\) es la función de distribución acumulativa de la distribución normal estándar y \(\beta\) son los parámetros del modelo, estimados a través de máxima verosimilitud.
• El modelo puede ser expresado de la siguiente manera: \(Y = X^{T}\beta+\epsilon\), donde \(\epsilon \sim N(0, 1)\).
• Las funciones logística y probit difieren en la manera como definen la función de distribución, mientras que la primera utiliza la función logística la segunda hace uso de la función de distribución acumulada de la normal estándar. Ambas funciones pueden ser comparadas en la siguiente figura:

Estadísticos descripvitos

library(tidyr)
datos %>% 
  gather(key = "variable", value = "valor", -r) %>% 
  group_by(variable, r) %>% 
  summarise(Promedio = round(mean(valor, na.rm = TRUE), digits = 2),
            `D. Estándar` = round(sd(valor, na.rm = TRUE), digits = 2),
            Mínimo = round(min(valor, na.rm = TRUE), digits = 2),
            Máximo = round(max(valor, na.rm = TRUE), digits = 2),
            Q1 = round(quantile(valor, prob = 0.25, na.rm = TRUE), digits = 2),
            Q2 = round(quantile(valor, prob = 0.5, na.rm = TRUE), digits = 2),
            Q3 = round(quantile(valor, prob = 0.75, na.rm = TRUE), digits = 2)) %>% 
  rename(Oxalato = r, Variable = variable) %>% 
  datatable()

## `summarise()` regrouping output by 'variable' (override with `.groups` argument)

Análisis Exploratorio

Datos ausentes

library(broom)
tidy(apply(datos, MARGIN = 2, is.na)) %>% 
  gather(key = "variable", value = "valor") %>% 
  mutate(valor = as.numeric(valor))  %>% 
  group_by(variable) %>% 
  summarise(Total_NAs = sum(valor))

## Warning: 'tidy.logical' is deprecated.
## See help("Deprecated")

## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

## `summarise()` ungrouping output (override with `.groups` argument)

## Warning: `...` is not empty.
## 
## We detected these problematic arguments:
## * `needs_dots`
## 
## These dots only exist to allow future extensions and should be empty.
## Did you misspecify an argument?

## # A tibble: 1 x 2
##   variable Total_NAs
##   <chr>        <dbl>
## 1 x                0

Densidades

library(ggplot2)
colores <- c("dodgerblue", "gray40")
datos %>% 
  rename(Oxalato = r) %>% 
  gather(key = "variable", value = "valor", -Oxalato) %>% 
  ggplot(data = ., aes(x = valor, fill = Oxalato)) +
  facet_wrap(~variable, scales = "free") +
  geom_density(alpha = 0.7) +
  scale_fill_manual(values = colores) +
  labs(x = "", y = "") +
  theme_light() +
  theme(legend.position = "bottom")

## Warning: Removed 2 rows containing non-finite values (stat_density).

Distribución condicional

• Las densidades condicionales son gráficos exploratorios que permiten dilucidar cómo es la probabilidad de “éxito” o “fracaso” respecto a variables numéricas. Es posible evidenciar en qué puntos (valores de x) se maximiza la probabilidad de “éxito”, que en este caso está ligado a la presencia (r = 1) de oxalatos en orina.

par(mfrow = c(2, 3))
cdplot(datos$r ~ datos$gravity, xlab = "Gravedad específica",
       ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$ph, xlab = "pH",
       ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$osmo, xlab = "Osmolaridad",
       ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$cond, xlab = "Conductividad",
       ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$urea, xlab = "Urea",
       ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$calc, xlab = "Calcio",
       ylab = "Oxalato", col = colores)

Boxplot comparativo

datos %>% 
  rename(Oxalato = r) %>% 
  gather(key = "variable", value = "valor", -Oxalato) %>% 
  ggplot(data = ., aes(x = Oxalato, y = valor, fill = Oxalato)) +
  facet_wrap(~variable, scales = "free") +
  geom_boxplot() +
  scale_fill_manual(values = colores) +
  labs(x = "", y = "") +
  theme_light() +
  theme(legend.position = "bottom")

## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

Estimación de los modelos

# Modelos Lineales Generalizados
mod_logit  <- glm(r ~ ., data = datos, family = binomial(link = "logit"))
mod_probit <- glm(r ~ ., data = datos, family = binomial(link = "probit"))

# -------- Validación LOOCV (Manual)
out_logi_0.5 <- NULL
out_prob_0.5 <- NULL
for(i in 1:nrow(datos)){
  out_logi_0.5[i] = predict(update(mod_logit, data = datos[-i, ]),
                   newdata = datos[i,], type = "response")
  out_prob_0.5[i] = predict(update(mod_probit, data = datos[-i, ]),
                   newdata = datos[i,], type = "response")
}

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

# -------- Validación LOOCV (cv.glm())

## Función de coste con cutoff = 0.5
coste_0.5 <- function(r, pi = 0) mean(abs(r-pi)> 0.5)

## LOOCV
library(boot)
cv_error_logi_0.5 <- cv.glm(data = datos %>% filter(!is.na(osmo) & !is.na(cond)),
                        glmfit = mod_logit, cost = coste_0.5)
cv_error_prob_0.5 <- cv.glm(data = datos %>% filter(!is.na(osmo) & !is.na(cond)),
                        glmfit = mod_probit, cost = coste_0.5)

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

• Observaciones: la ejecución automática de LOOCV con la función cv.glm() requiere una función de coste para calcular el error. El objeto devuelto por la función del paquete boot almacena el error con el nombre delta; al restar 1 menos el error (delta) se obtendrá la precisión o accuracy del modelo, que es exactamente el mismo valor obtenido manualmente.

• Función de coste: en el código hay una función de coste o pérdida para el límite igual a 0.5. Esta función puede ser expresada de la siguiente manera: \(error = \sum |r_i - p_i| > 0.5\). Donde \(r_i\) es el i-ésimo valor real y \(p_i\) es el i-ésimo valor predicho. En el siguiente código es posible evidenciar que se obtienen los mismos resultados de forma manual y con la función cv.glm().

Finalidad

Los modelos dicotómicos también sirven para la predicción, y según cierta condición, brindarán una probabilidad (entre 0-1) sobre un evento que así se consulte a través del modelo estimado.

9. Series temporales.

El análisis de las series temporales juega un rol esencial tanto para conocer un determinado fenómeno como para pronosticar hacia el futuro.

Las series de tiempo son observaciones sobre un determinado fenómeno efectuadas en el tiempo, en lapsos ojala equivalente, o con intervalos de igual valor.

Ejemplos: exportaciones mensuales, ventas diarias de un producto, casos semanales de sida, temperatura promedio diaria, tasa anual de mortalidad, numero mensual de divorcios.

Una serie de tiempo posee dos características esenciales:

Los valores están ordenados o presentados en el tiempo.

Existe una dependencia o correlación entre los valores de la serie en el tiempo.

Un análisis por series temporales podría buscar:

Descripción de la serie Adecuar un modelo o técnica estocástica Pronostico en un periodo h de la serie.

El análisis de la serie debe preguntarse sobre el tipo de serie que se está analizando, los tipos de datos, y el período de referencia en la adecuación del mejor modelo estocástico.

De igual forma, dependiendo de la serie, el pronóstico debe considerar ciertas restricciones.

El análisis de una serie de tiempo de proceder, sin ser restrictivo, procede a groso modo de la siguiente forma:

• Descripción de la serie temporal.
• Selección de datos: rango de modelización, rango de estimación, pronóstico.
• Estimación de los modelos.
• Medidas de rendimiento (set de entrenamiento).
• Estimación al set de validación.
• Medidas de rendimiento (set de validación).
• Selección del mejor método de estimación
• Pronóstico.

Predicción vs. Pronóstico

Para faciliar las cosas, hablemos de la regresión en un contexto transversal y otro longitudinal.

En un contexto transversal, nuestro puntos o sujetos de estudio son las personas en un momento ti. En un contexto longitudinal, nuestro puntos o sujetos de estudio son los peridos de todo un rango de períodos t.

Mientras que la predicción busca saber el efecto de los coeficientes y realizar una intrapolación, en el pronóstico buscamos lo opuesto: extrapolamos y no le damos importancia al valor de los coeficientes.

Ejemplo Pronóstico del nivel del CO2

Primero cargamos los datos y los transformados en ts (time series).

co2 <- read.csv("C:/Users/oscar/Desktop/R --- SAF/Tema 7/co2.csv", sep=";")

co2ts<-ts(co2$CO2, start = c(1959,1), frequency = 12)

La estructura de una serie ts:

print(co2ts)

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1959 333.13 332.09 331.10 329.14 327.36 327.29 328.23 329.55 330.62 331.40
## 1960 333.92 333.43 331.85 330.01 328.51 328.41 329.25 330.97 331.60 332.60
## 1961 334.68 334.17 332.96 330.80 328.98 328.57 330.20 331.58 332.67 333.17
## 1962 336.82 336.12 334.81 332.56 331.30 331.22 332.37 333.49 334.71 335.23
## 1963 337.95 338.00 336.37 334.47 332.46 332.29 333.76 334.80 336.00 336.63
## 1964 339.05 339.27 337.64 335.68 333.77 334.09 335.29 336.76 337.77 338.26
## 1965 341.47 341.31 339.41 337.74 336.07 336.07 337.22 338.38 339.32 340.41
## 1966 343.02 342.54 340.88 338.75 337.05 337.13 338.45 339.85 340.90 341.70
## 1967 344.28 343.42 342.02 339.97 337.84 338.00 339.20 340.63 341.41 342.68
## 1968 345.92 345.40 344.16 342.11 340.11 340.15 341.38 343.02 343.87 344.59
## 1969 347.38 346.78 344.96 342.71 340.86 341.13 342.84 344.32 344.88 345.62
## 1970 348.53 347.87 346.00 343.86 342.55 342.57 344.11 345.49 346.04 346.70
## 1971 349.93 349.26 347.44 345.55 344.21 343.67 345.09 346.27 347.33 347.82
## 1972 351.71 350.94 349.10 346.77 345.73                                   
##         Nov    Dec
## 1959 331.87 333.18
## 1960 333.57 334.72
## 1961 334.86 336.07
## 1962 336.54 337.79
## 1963 337.93 338.95
## 1964 340.10 340.88
## 1965 341.69 342.51
## 1966 342.70 343.65
## 1967 343.04 345.27
## 1968 345.11 347.07
## 1969 347.23 347.62
## 1970 347.38 349.38
## 1971 349.29 350.91
## 1972

Visualización de la serie

Veamos la serie como tal:

autoplot(co2ts, ts.colour = "blue", ts.linetype = "dashed")

¿Qué podemos decir? ¿Hay autocorrelación y estacionalidad en el tiempo?

Un gráfico de los autocorrelaciones totales nos puede ayudar:

autoplot(acf(co2ts, plot = FALSE))

## Warning: `mutate_()` is deprecated as of dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Vemos un estudio de los componentes de una serie:

autoplot(stl(co2ts, s.window = "periodic"), ts.colour = "blue")

Estimación de modelos de series de tiempo: ARIMA

Vamos a estimar un total de 6 modelos ARIMA y luego comprar al mejor entre ellos:

• ARIMA(0,1,2)(0,1,1)
• ARIMA(1,1,0)(2,1,0)
• ARIMA(1,1,2)(2,1,1)
• ARIMA(1,1,1)(2,1,1)
• ARIMA(1,1,2)(1,1,1)
• ARIMA(0,1,1)(0,1,1)
• ARIMA(1,1,0)(1,1,0)

arima1<- Arima(co2ts, order=c(0,1,2), seasonal=list(order=c(0,1,1),period=12))
arima2<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(2,1,0),period=12))
arima3<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(2,1,1),period=12))
arima4<- Arima(co2ts, order=c(1,1,1), seasonal=list(order=c(2,1,1),period=12))
arima5<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(1,1,1),period=12))
arima6<- Arima(co2ts, order=c(0,1,1), seasonal=list(order=c(0,1,1),period=12))
arima7<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(1,1,0),period=12))

Veamos las salidas del ARIMA(0,1,2)(0,1,1)

arima1

## Series: co2ts 
## ARIMA(0,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     ma2     sma1
##       -0.5126  0.0272  -0.7499
## s.e.   0.0813  0.0854   0.0920
## 
## sigma^2 estimated as 0.1041:  log likelihood=-38.86
## AIC=85.73   AICc=86.01   BIC=97.71

¿Qué nos dice la salida? ¿Cómo interpretamos los coeficientes?

Veamos cuál es el mejor modelo entre los 6 estimados. Utilizaremos el criterio de AIC y BIC:

AIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)

##        df       AIC
## arima1  4  85.72537
## arima2  4  93.70131
## arima3  7  90.55929
## arima4  6  89.02326
## arima5  6  88.62985
## arima6  3  83.82626
## arima7  3 106.19086

BIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)

##        df       BIC
## arima1  4  97.71422
## arima2  4 105.69016
## arima3  7 111.53978
## arima4  6 107.00654
## arima5  6 106.61312
## arima6  3  92.81789
## arima7  3 115.18250

Aquí se puede apreciar que los ajustes que mejor AIC y BIC presentan son aquellas que solo tienen componente de médias móviles y no tienen componente autorregresiva, siendo ARIMA(0,1,1)(0,1,1) el modelo que los test arrojan con un menor valor y por tanto con una mayor consideración.

Una vez estimados los modelos y elegido el mejor de ellos, en este caso ARIMA(0,1,1)(0,1,1), se procede a validarlo, esto es, ver la adecuidad del modelo seleccionado. Para ello en primer lugar se grafican los correlogramas de los residuos para comprobar que son ruido blanco:

ggtsdiag(arima6)

El modelo respecta los criterior de adecuación.

Pronóstico –> horizonte de proyección h

Finalizada la etapa de estimación, selección y validación de un modelo candidato, vamos a proyectar para los siguientes cuatro años siguientes, lo cuál son 48 meses, o un horizonte de proyección h=48.

forecast1<-forecast(arima6, level = c(95), h = 50)
autoplot(forecast1)

## Warning: `filter_()` is deprecated as of dplyr 0.7.0.
## Please use `filter()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Otras técnicas de estimación.

En ete caso solo se detalló una estimación por un modelo ARIMA, pero para series univariadas se podría optar por:

• Regresión en el tiempo
• Descomposición de una TS
• ETS
• Regresiones dinámicas
• Redes neuronales
• Entre otros…

Fin